%matplotlib widget
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Sampling from 2D distributions

Sampling uniformly from unit disc

To generate samples in the unit disc it is convenient to generate them in polar coordinate form. Unfortunately uniformly generating $r$ and $\theta$ does not do it:

def plot_circle(r, c='k'):
    plt.plot(r * np.cos(np.linspace(0,2*np.pi,121)), r * np.sin(np.linspace(0,2*np.pi,121)), c)

rs = np.random.random_sample(10000)
thetas = 2 * np.pi * np.random.random_sample(10000)
plt.figure()
plt.plot(rs * np.cos(thetas), rs * np.sin(thetas), '.', ms=1)
plt.axis('equal')
plot_circle(1)

Figure

These samples are clumped towards the middle because the radii are uniformly distributed, so there are roughly equal numbers of points with $0 < r < 0.1$ and $0.9 < r < 1$. But the areas on the plane corresponding to these two intervals in $r$ are very different.

plot_circle(0.1)
plot_circle(0.9)

The area occupied by the interval between $r$ and $r+dr$ is proportional to $r$, so to generate points uniformly we want $p(r) \propto r$. The integral of this is $\frac{1}{2}r^2$, so the normalized pdf is $p(r) = 2r$ with the cdf $P(r) = r^2$. This means if we have a random number $\xi$ that's uniform on the unit interval, we can compute $r = \sqrt{\xi}$ and then $r \sim p(r)$.

rs = np.sqrt(np.random.random_sample(10000))
plt.figure()
plt.plot(rs * np.cos(thetas), rs * np.sin(thetas), '.', ms=1)
plt.axis('equal')
plot_circle(1)

Figure

Hooray! Now the points are uniformly distributed.

Generating points on the unit hemisphere

In rendering, the random selections we are making are often directions, or point on the sphere, or hemisphere if we are talking about directions on one side of some surface. The most basic distribution here is the uniform distribution. How can we generate points that are uniformly distributed on the hemisphere? As with the disc, uniformly choosing $\theta$ and $\phi$ is not going to work.

def plot_hemisphere_points(xs, ys, zs):
    fig = plt.figure(figsize=(5,5))
    ax = fig.add_subplot(111, projection='3d')
    ax.set_aspect('equal')
    plot_circle(1)
    ax.plot(xs, ys, zs, '.', ms=1)
    ax.set_xlim(-1,1)
    ax.set_ylim(-1,1)
    ax.set_zlim(0,2)

thetas = np.pi/2 * np.random.random_sample(4000)
phis = 2 * np.pi * np.random.random_sample(4000)
xs = np.sin(thetas) * np.cos(phis)
ys = np.sin(thetas) * np.sin(phis)
zs = np.cos(thetas)
plot_hemisphere_points(xs, ys, zs)

Figure

The situation is quite similar, but the details are different: the surface area of a strip from $\theta$ to $\theta + d\theta$ is proportional to $\sin\theta$. The cdf is $1 - \cos\theta$, so it is already normalized. Thus we can use the transformation $\theta(\xi) = \cos^{-1}(1 - \xi)$ and we have $\theta \sim \sin\theta$.

xis = np.random.random_sample(4000)
thetas = np.arccos(1 - xis)
phis = 2 * np.pi * np.random.random_sample(4000)
xs = np.sin(thetas) * np.cos(phis)
ys = np.sin(thetas) * np.sin(phis)
zs = np.cos(thetas)
plot_hemisphere_points(xs, ys, zs)

Figure

If you look at how we are computing the $z$ coordinate here, we are computing an arccosine followed by a cosine, which is a little silly. We could instead just set $z = 1 - \xi$ and then compute $\sin\theta$ as $\sqrt{1 - z^2}$.

xis = np.random.random_sample(4000)
zs = 1 - xis
phis = 2 * np.pi * np.random.random_sample(4000)
xs = np.sqrt(1 - zs**2) * np.cos(phis)
ys = np.sqrt(1 - zs**2) * np.sin(phis)
plot_hemisphere_points(xs, ys, zs)

Figure

One last very fundamental distribution on the hemisphere is the cosine distribution, in which we'd like the density on the hemisphere to be $p(\omega) = \cos\theta / \pi$. (Here $\theta$ depends on $\omega$.) The area of a $d\theta$ slice is still proportional to $\sin\theta$, so if the points are cosine distributed on the hemisphere we will find $p(\theta) \propto \cos\theta \sin\theta$. Let's work this out out interactively in class....

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

Hopefully we came to the conclusion that the proper cdf is $P(\theta) = \frac{1}{2}(1 - \cos 2\theta)$; solving, $\theta = \frac{1}{2} \cos^{-1}(1 - 2\xi)$.

xis = np.random.random_sample(4000)
thetas = 1/2 * np.arccos(1 - 2*xis)
phis = 2 * np.pi * np.random.random_sample(4000)
xs = np.sin(thetas) * np.cos(phis)
ys = np.sin(thetas) * np.sin(phis)
zs = np.cos(thetas)
plot_hemisphere_points(xs, ys, zs)

Figure

An interesting feature of this point set is that if you project it down onto the unit disc (for instance, by viewing the plot from straight overhead), it looks uniform. This is because an area $dA$ on the sphere projects to an area $dA \cos\theta$ on the equatorial plane, so the cosine-weighted distribution is the correct one.

This observation leads to a question: can we simplify the formulas we just derived to something that looks more like sample-disc-then-unproject? Well, again the observation of code that computes a trig function of something it just computed using an inverse trig function suggests we can do better. Our sampling rule could be written as: $$ \begin{aligned} 1 - 2\xi &= \cos 2\theta = 1 - 2 \sin^2 \theta \\ \xi &= \sin^2\theta \\ 1-\xi &= \cos^2\theta \end{aligned} $$ ...so we can compute $\sin\theta$ (the distance from the $z$ axis) as $\sqrt\xi$, and then $\cos\theta$, the height above the $x-y$ plane, is $\sqrt{1-\xi}$.

xis = np.random.random_sample(4000)
rs = np.sqrt(xis)
phis = 2 * np.pi * np.random.random_sample(4000)
xs = rs * np.cos(phis)
ys = rs * np.sin(phis)
zs = np.sqrt(1 - xis)
plot_hemisphere_points(xs, ys, zs)

Figure

Note that now, the first 5 lines are exactly the same as the code to sample the unit disc.